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We obtain exact results on interaction quenches in the ID Bose gas described by the integrable 
Lieb-Liniger model. We show that in the long time limit integrability leads to significant deviations 
from the predictions of the grand canonical ensemble and a description within the generalized Gibbs 
ensemble (GGE) is needed. For a non-interacting initial state and arbitrary final interactions, we 
find that the presence of infinitely many conserved charges generates a non-analytic behavior in the 
equilibrated density of quasimomenta. This manifests itself in a dynamically generated Friedel-like 
oscillation of the non-local correlation functions with interaction dependent oscillation momenta. 
We also exactly evaluate local correlations and the generalized chemical potentials within GGE. 

PACS numbers: 



Whether and how an isolated quantum system equili- 
brates or thermalizes are fundamental questions in un- 
derstanding non-equilibrium dynamics. The answers can 
also shed light on the applicability of quantum statisti- 
cal mechanics to closed systems. While these questions 
are very hard to study experimentally in the condensed 
matter setup, they have become accessible in ultracold 
quantum gases due to the recent experimental advances 
[1]. Thanks to their unprecendented tunability, ultracold 
atomic systems allow for the study of non-equilibrium 
quantum dynamics of almost perfectly isolated strongly 
correlated many-body systems in a controlled way. These 
experiments [2, 3] triggered a revival of theoretical studies 
on issues of thermalization [4] . Fundamental questions in- 
clude whether stationary values of local correlation func- 
tions are reached in a system brought out of equilibrium, 
and if so, how they can be characterized. Can conven- 
tional statistical ensembles describe the state? Is there 
any kind of universality in the steady state and the way 
it is approached? 

The absence of thermalization of a ID bosonic gas re- 
ported in [3] brought to light the special role of inte- 
grability. The lack of thermalization was attributed to 
the fact that the system was very close to an integrable 
one, the Lieb-Liniger (LL) model [5], which is the sub- 
ject of our Letter. The dynamics of integrable systems is 
highly constrained by the presence of a large number of 
so-called conserved charges in addition to the total par- 
ticle number, momentum, and energy. Thus integrable 
systems are not expected to thermalize. Following the 
ideas of the subjective statistical mechanics of Jaynes [6], 
the so-called generalized Gibbs ensemble (GGE) was pro- 
posed [7] to capture the long-time behavior of integrable 
systems brought out of equilibrium. This ensemble is the 
least biased statistical representation of the system once 
the conserved charges {Qm} are taken into account. The 
density matrix is 



where the generalized "chemical" potentials {/3m} are 
fixed by the expectation values (Qm), and Zqq^ = 



Tr 



The GGE was tested and its drawbacks were analyzed 
by various numerical and analytical approaches [8] . How- 
ever, with a few exceptions [9, 10], only models that 
can be mapped to quadratic bosonic or fermionic sys- 
tems have been considered, where the conserved charges 
are given by the mode occupation numbers. While some 
of these models are paradigmatic, like the Ising or Lut- 
tinger models, a prominent class of non-trivial integrable 
systems has not been sufficiently explored, namely those 
solvable by the Bethe Ansatz. In these models, the con- 
served charges are local and cannot be expressed as mode 
occupations. 

In this Letter, we derive experimentally testable pre- 
dictions for the long time behavior of the system after 
an interaction quench [11, 12] in the LL model by com- 
bining its Bethe Ansatz solution and GGE. For a non- 
interacting initial state and arbitrary final interactions, 
the infinitely many conserved charges result in Friedel- 
like oscillations in the two-point functions, in striking 
contrast to the grand canonical ensemble (GCE). 

The model. — The LL model describes a system of 
identical bosons in ID interacting via a Dirac-delta po- 
tential. The Hamiltonian in second quantized formulation 
is given by [5] 



H = 



(^dj^dj + ci>^^i;i>'^ , (2) 



Pgge 



Zgge 



(1) 



where c > in the repulsive regime we wish to study, 
and for brevity we have set H = 1 and the boson mass 
to be equal to 1/2. The dimensionless coupling constant 
is given by 7 = c/n, where n = N/L is the density of 
the gas. In cold atom experiments 7 is a function of the 
3D scattering length and the ID confinement [3, 13]. The 
exact spectrum and thermodynamics of the model can be 
obtained via Bethe Ansatz [5, 14]. Each eigenstate of the 
system with N particles on a ring of circumference L is 
characterized by a distinct set of quantum numbers {Ij} 



2 



that are integers (half-integers) for A'' odd (even). The 
wave function can be expressed in terms of N quasimo- 
menta {Xj} that satisfy a set of algebraic equations 



N 
k=l 



Afe) = 2ttL 



(3) 



where 9{X) = 2 arctan(A/c). The expectation values of 
the local conserved charges can be computed as 



Qni — {Qr> 



(4) 



in particular, the energy of a given state is simply E ~ 
Q2 — The wave function is identically zero if 

any two of the {Ij} coincide, which is reminiscent of the 
Pauli principle for fcrmions. In the thermodynamic limit, 
a mixed state corresponding to thermal equilibrium [15] 
is captured by a filling fraction < // < 1 in the space 
of quantum numbers. It has been shown in Ref. [10] that 
GGE also implies that such a description is possible. For 
calculations, it is more convenient to define a function 
/(A) in terms of the quasimomenta. Unlike the case of 
free fermions, where the relation between // and /(A) is 
trivial, here all quasimomenta are coupled to each other 
by Eq. (3). Thus the density of occupied quasimomenta 
p^P^(A) is not independent of /(A) but satisfies the inte- 
gral equation [5] 



P(A) 



2^ 27 

p(p)(A)=/(A)p(A) 



(5a) 
(5b) 



with the kernel ip{X) = 2c/(A^ + c^). Here p{X) is the 
maximal allowed density of quasimomenta and p^P' (A) is 
normalized as / dA p^^^ (A) = n. 

Conserved charges and the problem of moments. — 
The simplest way to bring a system out of equilibrium 
is a sudden change of one of its parameters, a quantum 
quench. In a cold atom setting such a quench could be 
achieved by a rapid change of the transverse confinement 
or the scattering length. We will compute the predictions 
of the GGE for a sudden quench of the interaction param- 
eter, and compare them to those of the GCE. In order to 
describe the final state in terms of p*-P-'(A), one needs to 
find the expectation values of the conserved charges Qm 
right after the quench. For c > 0, all solutions of Eq. (3) 
are real [14], thus finding the density p^^\X) which repro- 
duces these charges is equivalent to solving the problem 
of moments defined by 



L I dAp(P'(A)A 



(6) 



The first few conserved charge operators Qm can be 
written in terms of the field operator ijj{x) as Qo = 
Jdx^^^, J Qi — —i Jdxtp^dxip, and Q2 = H is the 



Hamiltonian given by Eq. (2). Except for Qo and Qi, the 
charges depend on c and their expectation values change 
during the quench according to their c-dependence [11]. 
Unfortunately, the corresponding expressions are not 
known explicitly for the operators Qm for to > 4 [16]. 

In the thermodynamic limit, the even conserved 
charges are expected to have the form 



Q2 



dx 



{deriv.} + Am c" 



m+l 



(7) 



where {deriv.} denotes unknown terms which involve 
products of spatial derivatives of tp{x) and tp^{x) and 
: : denotes normal ordering. The constants Am can be 
calculated comparing Eq. (7) to the semiclassical limit 
(7^1) where the expressions for the conserved charges 
are known [17]. The constants Am satisfy a certain re- 
cursion relation, which can be solved to obtain Am ~ 
2™(2to - 1)!!/(to + 1)! [23]. Note that assuming an even 
p^P)(A), as in the ground state, we keep only the even 
charges. 

Universal semicircle solution. — If the initial state is 
in a pure non-interacting BEG (although we expect our 
results to be also valid for small initial interactions) , then 
the expectation values of the unknown derivative terms 
in Eq. (7) vanish. Thus the conserved charges are given 
by Q2ra/L « A„c"^(: (^/'t^)™+i :) « A^n^^+i-y™, where 
in the last step we used the fact that the local correlators 
(l^tm^m^ — >■ as 7 — 0. We are left with the following 
problem of moments: 



dAp(p)(A) A" = 2' 



(m-l- 1)! 



7" , (8) 



which has the unique solution 



p(p)(A) 



on the interval [—A*, A*] and zero otherwise (see inset of 
Fig. 1). Using Eqs. (5) we can determine /(A): 



/(A) 



(1 - {X/X^f) 



f + y(V-(A/A*)')'+7 



^ - (A/A. 



(10) 

We emphasize that the existence of the sharp edge of sup- 
port in Eq. (9) is a direct consequence of our accounting 
for an infinite number of conserved charges. Keeping any 
finite number within GGE would smear the sharp edge. 

Correlation functions in the final state. — Knowing 
the filling fraction /(A) allows us to calculate some of 
the correlation functions. First we calculate local corre- 
lators exactly using the results of Ref. [18] which give 
analytic expressions for the local two and three-point 
correlators for arbitrary states that are captured by a 
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FIG. 1: (Color online) Quench from a non-interacting ini- 
tial state to arbitrary final interactions. Local correlations 
g-i and as functions of the coupling 7, calculated from 
the generalized Gibbs ensemble (solid, black) and from the 
grand canonical ensemble (dots, red). Inset: density of quasi- 
momenta p'P'(A) in the two ensembles for 7 = 100. 




FIG. 2: (Color online) Equal time Green's function for quasi- 
momenta distributed with the semi-circle function (dots) and 
with the thermal distribution (solid, red) when 7 = 400. In- 
set: dynamically generated Friedel-like oscillation momentum 
k, as a function of 7 in units of fcp = mi (solid, black) with 
the large 7 asymptotic result (dashed, blue). 



continuous /(A). Wc compute 52 = (: (i/'W-')^ ■)/"-^ a-nd 
53 = (: (?/>t?/,)3 :) /n^ both for both for GGE and the GCE 
by using the appropriate /(A). In the latter only the en- 
ergy and the particle densities are fixed to be the same 
as for the GGE, which determines the temperature and 
the chemical potential. The results are shown in Fig. 1. 
For interactions not too strong, the GGE results are very 
close to the values for GCE which slightly underestimates 
52 and overestimates g^. The deviation is bigger for 
and grows with 7. Both GGE and GCE results are con- 
sistent with the outcomes of recent numerical simulation 
[11]. We conclude that the local correlators above are 
too simple to effectively distinguish between the various 
ensembles. 

Note, however, that the densities of filled states /9^p^(A) 
are very different for the GGE and the GCE, see the in- 
set of Fig. 1. Remarkably, the semicircle distribution ex- 
hibits a non-analytic square root singularity at the edges, 
which has immediate consequences for the behavior of 
the nonlocal correlation functions. In analogy with the 
Friedel oscillations, which appear due to the existence of 
a sharp Fermi surface, we expect that the non-analyticity 
of yO^P^(A) at A* will result in the oscillating behavior of 
nonlocal correlation functions. Unlike the period of con- 
ventional Friedel oscillations which depends only on the 
average density, the period of these dynamically gener- 
ated Friedel oscillations also depends on the interaction 
strength. The characteristic momentum fc* which plays 
the role of the Fermi momentum can be calculated using 
standard techniques [5, 14], and is plotted in the inset 
of Fig. 2. While we cannot explicitly calculate the nonlo- 
cal correlation functions for arbitrary final interactions, 
below we present detailed results for strong final interac- 
tions which support this picture. 



Strongly interacting final state. — For large coupling 
the system is effectively in the fermionized regime, since 
the largest quasimomentum A* scales only as ^/c and 
thus the convolution term in Eq. (5a) depending on the 
kernel (p{X) can be neglected. The function /(A) takes 
the form 



7^00 ^ y 



(11) 



Bosonic correlation functions can be calculated by first 
fermionizing the field operators using Jordan- Wigner 
strings, ip{x) = exp[iTT J'^^ipl{z)'4>F{z)dz]'tljF{x), and 

then exploiting free fermionic correlators of ipp in the 
large 7 regime. Let us consider the equal time Green's 
function 



G(x) - (V3^(a;)7^(0))sc 



(12) 



where SC denotes averaging over the semi-circle distri- 
bution of Eq. (11). We proceed by introducing a lattice 
discretization a <^ l/n,x — (m + l)a and recasting the 
correlator as 

G(x) = Ul{x) Y[ [l-2 4,l{ka)Mka)\ ^fW) • 

\ k<m I 

(13) 

The long chain of operators is now amenable to a 
Wick expansion using as a building block the fermionic 
two point function given by the Fourier transform 
of the distribution (11), Gff{x) = J ^fWe'^'^ = 
2n Ji(A*a;)/(A*a;), where Ji{x) is the Bessel function 
of the first kind. The Wick expansion of Eq. (13) can 
be recast as a Fredholm-like determinant since the cor- 
relator of 2n fields can be represented as a deter- 
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FIG. 3: (Color online) Dimensionless inverse temperature 
n^P2 (dots) as a function of 7, shown together with the grand 
canonical prediction (solid, red) and analytical 7 3> 1 re- 
sult « 1/(87) (dashed, blue). Inset: generalized dimensionless 
chemical potential /?2mAj"' as a function of m for 7 = 1 with 
large m result (15) indicated. 



minant of size n, and we need to sum over all dis- 
cretizations of the interval (0, x) which can be accom- 
pUshed via a determinant expansion [19]. Finally, we 
have G{x) = lima_>.o detm-i-i(2aB — A)/(2a) with Aij = 
diag{0, 1, . . . , 1} and Bij = Gff {{i — j)a) for i,j > 1. 
The limit a — > is taken numerically. For the semicir- 
cle distribution we find that the correlation function can 
be well fitted using the Bessel function with an exponen- 
tially decaying term, 2n Ji(A*a;)/(A*a;)e~"'^. The value 
of a changes from 2.25 to 2 as 7 is varied from 100 to 
infinity. A field theoretical understanding of the behavior 
of correlation functions, possibly in the spirit of models 
going beyond the Luttinger Liquid in ID [20] remains to 
be achieved. 

We can also compute the density-density correlation 
function (72(2^) ~ {x)tp'^ {0)t[}{0)rlj{x)) /n^ for large final 
7 using the first few terms of the infinite series given 
in Ref. [21]. In the large 7 limit the leading order for 

arbitrary /(A) is given by 32 (x) « 1 - (/ f^/(A)e''^-)' . 
For the GGE case we use the /(A) given in Eq. (11). 
The result, 32(2;) ~ 1 — 2nJ^{X^x)/{X^x)^ , is osciUatory 
and decays as a power law ^ l/.i'^, in contrast to the 
exponential decay in the thermal state. 

GGE generalized chemical potentials. — The GGE 
density matrix (1) minimizes the generalized free en- 
ergy, J"/(fcBT) = X]m/5™Q™ ~ "^/^B- Expressing {Qm} 
and the entropy S in terms of p(A) and /(A), we can 
obtain a nonlinear integral equation for the filling frac- 
tion /(A) = l/(e'^''^^ + 1) following the standard proce- 
dure [10, 15, 22]: 

£(A)+ / ^^(A-A')log(l + e-^(^')) =^/3™A". 

rn 

(14) 



Here /32 = /? = l/{k'QT) is the inverse temperature. Note 
that keeping only the particle number Qq = N and en- 
ergy Q2 = E, one recovers the usual GCE in Eq. (1) and 
the standard Yang- Yang equation [15] in Eq. (14). 

Since we know £(A) = log[l//(A) — 1] from the so- 
lution of the problem of moments given by Eq. (10), 
the generalized chemical potentials {f32m} can be deter- 
mined by expanding the left hand side of Eq. (14) in pow- 
ers of A/ A*. The convolution term is a smooth bounded 
function, so for large m its expansion coefficients decay 
exponentially fast. However, the coefficients of e(A) de- 
cay only algebraically, since it diverges at A ^ ±A* as 
e(A) ~ — i log[l— (A/A*)^]. This divergence results in the 
large m behavior of the chemical potentials as 



" 2mA2'" 2m(4n27)™ ' ^ ' 

None of them diverges but since there are infinitely many 
of them, their collective behavior renders £(A) divergent 
for A > A*. This creates a "Fermi edge" by constraining 
p'^P^(A) to have a finite range of support. This mechanism 
is markedly different from the case of the ground state, 
where the divergence of a finite number of chemical po- 
tentials, as the temperature goes to zero, results in a step 
distribution function. 

Summary. — We studied the large time behavior of 
the post interaction quench state in the Lieb-Liniger 
model using analytic techniques by combining the gen- 
eralized Gibbs ensemble and Bethe Ansatz integrability 
of the model. For a non-interacting initial state and ar- 
bitrary final interactions, the infinitely many conserved 
charges result in Friedel-like oscillations in the two-point 
functions with interaction dependent momentum. This 
is in striking contrast to the behavior of the thermalized 
system and provides a smoking gun signal for cold atom 
experiments. We also exactly evaluated local correlations 
and the generalized chemical potentials within GGE. 
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